package fem2.material;

import static java.lang.Math.log;
import static java.lang.Math.pow;
import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import fem2.Element;
import fem2.NotImplementedError;
import fem2.element.StructuralElement;
import fem2.enu.State;

public class ModifiedStVenantKirchhoffMM extends ElasticMaterialModel {

	public ModifiedStVenantKirchhoffMM(double E, double nu, double thickness, double rho, State ss) {
		super(E, nu, thickness, rho, ss);
	}

	@Override
	public IMatrix getConstitutiveMatrix(Element e, int k) {
		int dim = e.getMeshPart().getDimension();
		int ntrs = dim * (dim + 1) / 2;
		double[] xi = e.getIntegrationPoint(k);
		IVector epsilon = ((StructuralElement) e).computeStrain(xi);
		IMatrix C = Array1DMatrix.createRowMajor(ntrs, ntrs);
		computeC(C, epsilon);
		return C;
	}

	@Override
	public void Initialize(Element e) {
	}

	@Override
	public void InitializeSolutionStep(Element e) {
		/**
		 * Calculate first stress for the element
		 */
		StructuralElement se = (StructuralElement) e;
		for (int i = 0; i < se.getNumIntegrationPoints(); i++) {
			IVector sigma = se.getStress(i);
			double[] xi = se.getIntegrationPoint(i);
			IVector epsilon = se.computeStrain(xi);
			computeStress(sigma, epsilon);
		}
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
		/*
		 * Calculate material response
		 */
		StructuralElement se = (StructuralElement) e;
		for (int i = 0; i < se.getNumIntegrationPoints(); i++) {
			IVector sigma = se.getStress(i);
			double[] xi = se.getIntegrationPoint(i);
			IVector epsilon = se.computeStrain(xi);
			computeStress(sigma, epsilon);
		}
	}

	@Override
	public void FinalizeSolutionStep(Element e) {
	}

	@Override
	public void CalculateOutputData(Element e) {
	}

	public void computeStress(IVector s, IVector e) {
		switch (state) {
		case PLANE_STRAIN:
			double e11 = e.get(0);
			double e22 = e.get(1);
			double e12 = e.get(2);

			double s11 = ((kappa + 2 * e22 * kappa)
					* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + 4 * e11 * mu
					* (1 + 2 * e22 + e11 * (2 + 4 * e22) - pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1);
			double s22 = ((kappa + 2 * e11 * kappa)
					* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + 4 * e22 * mu
					* (1 + 2 * e22 + e11 * (2 + 4 * e22) - pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1);
			double s12 = -(e12
					* (kappa * log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + 2 * mu
							* (-1 - 2 * e22 - 2 * e11 * (1 + 2 * e22) + pow(e12, 2))) * pow(2 + 4
					* e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1));
			s.set(0, s11);
			s.set(1, s22);
			s.set(2, s12);
			return;
		case PLANE_STRESS:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		case THREE_DIMENSIONAL:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		default:
			throw new Error("invalid state");
		}
	}

	public void computeC(IMatrix c, IVector e) {
		switch (state) {
		case PLANE_STRAIN:
			double e11 = e.get(0);
			double e22 = e.get(1);
			double e12 = e.get(2);

			double c11 = -((4 + 8 * e22)
					* ((kappa + 2 * e22 * kappa)
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + 4 * e11
							* mu * (1 + 2 * e22 + e11 * (2 + 4 * e22) - pow(e12, 2))) * pow(2 + 4
					* e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -2))
					+ pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* (4 * e11 * (2 + 4 * e22) * mu + 4 * mu
							* (1 + 2 * e22 + e11 * (2 + 4 * e22) - pow(e12, 2)) + (2 + 4 * e22)
							* (kappa + 2 * e22 * kappa)
							* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1));
			double c12 = -((4 + 8 * e11)
					* ((kappa + 2 * e22 * kappa)
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + 4 * e11
							* mu * (1 + 2 * e22 + e11 * (2 + 4 * e22) - pow(e12, 2))) * pow(2 + 4
					* e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -2))
					+ pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* (4 * e11 * (2 + 4 * e11) * mu + 2 * kappa
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + (2 + 4 * e11)
							* (kappa + 2 * e22 * kappa)
							* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1));
			double c13 = 4
					* e12
					* ((kappa + 2 * e22 * kappa)
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + 4 * e11
							* mu * (1 + 2 * e22 + e11 * (2 + 4 * e22) - pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -2)
					+ pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* (-8 * e11 * e12 * mu - 2 * e12 * (kappa + 2 * e22 * kappa)
							* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1));
			double c22 = -((4 + 8 * e11)
					* ((kappa + 2 * e11 * kappa)
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + 4 * e22
							* mu * (1 + 2 * e22 + e11 * (2 + 4 * e22) - pow(e12, 2))) * pow(2 + 4
					* e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -2))
					+ pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* (4 * (2 + 4 * e11) * e22 * mu + 4 * mu
							* (1 + 2 * e22 + e11 * (2 + 4 * e22) - pow(e12, 2)) + (2 + 4 * e11)
							* (kappa + 2 * e11 * kappa)
							* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1));
			double c23 = 4
					* e12
					* ((kappa + 2 * e11 * kappa)
							* log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + 4 * e22
							* mu * (1 + 2 * e22 + e11 * (2 + 4 * e22) - pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -2)
					+ pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* (-8 * e12 * e22 * mu - 2 * e12 * (kappa + 2 * e11 * kappa)
							* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1));
			double c33 = -4
					* pow(e12, 2)
					* (kappa * log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + 2 * mu
							* (-1 - 2 * e22 - 2 * e11 * (1 + 2 * e22) + pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -2)
					- (kappa * log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)) + 2 * mu
							* (-1 - 2 * e22 - 2 * e11 * (1 + 2 * e22) + pow(e12, 2)))
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					- e12
					* pow(2 + 4 * e22 + e11 * (4 + 8 * e22) - 2 * pow(e12, 2), -1)
					* (4 * e12 * mu - 2 * e12 * kappa
							* pow(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2), -1));

			c.set(0, 0, c11);
			c.set(0, 1, c12);
			c.set(0, 2, c13);
			c.set(1, 0, c12);
			c.set(1, 1, c22);
			c.set(1, 2, c23);
			c.set(2, 0, c13);
			c.set(2, 1, c23);
			c.set(2, 2, c33);
			break;
		case PLANE_STRESS:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		case THREE_DIMENSIONAL:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		default:
			throw new Error("invalid state");
		}
	}

	@Override
	public double computeEnergyAt(IVector e) {
		switch (state) {
		case PLANE_STRAIN:
			double e11 = e.get(0);
			double e22 = e.get(1);
			double e12 = e.get(2);

			/*
			 * TODO: check the equation below
			 */
			return mu * pow(e11, 2) + (mu * pow(e12, 2)) / 2. + mu * pow(e22, 2)
					+ (kappa * pow(log(1 + 2 * e11 + 2 * e22 + 4 * e11 * e22 - pow(e12, 2)), 2))
					/ 8.;
		case PLANE_STRESS:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		case THREE_DIMENSIONAL:
			/*
			 * TODO
			 */
			throw new NotImplementedError();
		default:
			throw new Error("invalid state");
		}
	}

	@Override
	public double computeEnergyAt(IVector e, IVector s) {
		/*
		 * TODO: check if there is better equation
		 */
		return computeEnergyAt(e);
	}
}
